--- permalink: /spectre/spatial-simple ---
Here we provide instructions for analysing spatial data after cell segmentation as part of Spectre’s simple segmentation and spatial analysis workflow. This workflow is split into 3x scripts, which can be run after performing cell segmentation:
If you go to https://github.com/ImmuneDynamics/Spectre you can download the repository, including the analysis workflow scripts.
You can find the simple spatial analysis workflow scripts in this folder.
You will need a folder called ‘data’ containing a folder called ‘ROIs’. Each ROI should be a folder with a unique name, and the TIFF files for that ROI should be in the corresponding folder.
Within the ‘data’ folder, you will need a folder called ‘masks’. Each of your mask types should be stored here – with the ROI name, and then some pattern that can separate which type of mask it is.
You will also need a ‘metadata’ folder containing a CSV file called ‘sample.metadata’. In this file you should include a column called ‘ROI’ with the name of each ROI, and then as many other columns as you like containing metadata for each ROI. We recommend including ‘Sample’ (e.g. if multiple ROIs are taken from one sample, you can note this here), ‘Group’ (the experiment group each ROI belongs to), and ‘Batch’ (the sample preparation/acquisition batch for the ROI – this could be some kind of reference, or a date).
…
###################################################################################
### Spectre: spatial 1 - add masks and extract cellular data
###################################################################################
First, load the Spectre package and associated packages.
### Load libraries
library('Spectre')
Spectre::package.check(type = 'spatial')
Spectre::package.load(type = 'spatial')
Next, set directories.
### Set PrimaryDirectory
dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
getwd()
PrimaryDirectory <- getwd()
PrimaryDirectory
### Set InputDirectory (ROI TIFFs)
setwd(PrimaryDirectory)
setwd("../data/ROIs/")
InputDirectory <- getwd()
InputDirectory
### Set MaskDirectory (ROI mask TIFFs)
setwd(PrimaryDirectory)
setwd("../data/masks")
MaskDirectory <- getwd()
MaskDirectory
### Create output directory
setwd(PrimaryDirectory)
dir.create("Output 1 - add masks")
setwd("Output 1 - add masks")
OutputDirectory <- getwd()
OutputDirectory
###################################################################################
### Check ROIs and TIFFs
###################################################################################
### Initialise the spatial data object with channel TIFF files
setwd(InputDirectory)
rois <- list.dirs(full.names = FALSE, recursive = FALSE)
as.matrix(rois)
## [,1]
## [1,] "ROI002"
## [2,] "ROI004"
## [3,] "ROI006"
## [4,] "ROI008"
## [5,] "ROI010"
## [6,] "ROI012"
### Check channel names
tiff.list <- list()
for(i in rois){
setwd(InputDirectory)
setwd(i)
tiff.list[[i]] <- list.files(getwd())
}
t(as.data.frame(tiff.list))
## [,1] [,2] [,3] [,4]
## ROI002 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI004 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI006 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI008 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI010 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI012 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## [,5] [,6] [,7] [,8]
## ROI002 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI004 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI006 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI008 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI010 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI012 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
###################################################################################
### Read in TIFF files and create spatial objects
###################################################################################
### Read in ROI channel TIFFs
setwd(InputDirectory)
spatial.dat <- read.spatial.files(dir = InputDirectory)
As the files are read in, you will see feedback like this:
## Reading TIFFs from:/Users/thomasa/OneDrive - The University of Sydney (Staff)/Library/Github (public)/Spectre/workflows/Spatial - simple/data/ROIs
## ...with extent correction
## ...with y-axis orientation flipping
## Reading ROI: ROI002
## -- reading single band in TIFF:CD11b_Sm149.tiff
## -- reading single band in TIFF:CD20_Dy161.tiff
## -- reading single band in TIFF:CD3_Er170.tiff
## -- reading single band in TIFF:CD4_Gd156.tiff
## -- reading single band in TIFF:CD45_Sm152.tiff
## -- reading single band in TIFF:CD8a_Dy162.tiff
## -- reading single band in TIFF:DNA1_Ir191.tiff
## -- reading single band in TIFF:DNA3_Ir193.tiff
### Check results
str(spatial.dat, 3)
## List of 6
## $ ROI002:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS : list()
## .. ..@ DATA : list()
## $ ROI004:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS : list()
## .. ..@ DATA : list()
## $ ROI006:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS : list()
## .. ..@ DATA : list()
## $ ROI008:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS : list()
## .. ..@ DATA : list()
## $ ROI010:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS : list()
## .. ..@ DATA : list()
## $ ROI012:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS : list()
## .. ..@ DATA : list()
spatial.dat[[1]]@RASTERS
## class : RasterStack
## dimensions : 501, 500, 250500, 8 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 0, 500, 0, 501 (xmin, xmax, ymin, ymax)
## crs : NA
## names : CD11b_Sm149, CD20_Dy161, CD3_Er170, CD4_Gd156, CD45_Sm152, CD8a_Dy162, DNA1_Ir191, DNA3_Ir193
## min values : 0, 0, 0, 0, 0, 0, 0, 0
## max values : 33, 779, 41, 40, 734, 109, 1670, 3007
###################################################################################
### Read in masks files
###################################################################################
### Define cell mask extension for different mask types
setwd(MaskDirectory)
all.masks <- list.files(pattern = '.tif')
as.matrix(all.masks)
## [,1]
## [1,] "ROI002_Cell_mask.tiff"
## [2,] "ROI002_Cytoplasm_mask.tiff"
## [3,] "ROI002_Nuclei_mask.tiff"
## [4,] "ROI004_Cell_mask.tiff"
## [5,] "ROI004_Cytoplasm_mask.tiff"
## [6,] "ROI004_Nuclei_mask.tiff"
## [7,] "ROI006_Cell_mask.tiff"
## [8,] "ROI006_Cytoplasm_mask.tiff"
## [9,] "ROI006_Nuclei_mask.tiff"
## [10,] "ROI008_Cell_mask.tiff"
## [11,] "ROI008_Cytoplasm_mask.tiff"
## [12,] "ROI008_Nuclei_mask.tiff"
## [13,] "ROI010_Cell_mask.tiff"
## [14,] "ROI010_Cytoplasm_mask.tiff"
## [15,] "ROI010_Nuclei_mask.tiff"
## [16,] "ROI012_Cell_mask.tiff"
## [17,] "ROI012_Cytoplasm_mask.tiff"
## [18,] "ROI012_Nuclei_mask.tiff"
mask.types <- list('cell.mask' = '_Cell_mask.tiff')
mask.types
## $cell.mask
## [1] "_Cell_mask.tiff"
### Read in masks
for(i in names(mask.types)){
spatial.dat <- do.add.masks(dat = spatial.dat,
mask.dir = MaskDirectory,
mask.pattern = mask.types[[i]],
mask.label = i)
}
str(spatial.dat, 3)
## List of 6
## $ ROI002:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI004:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI006:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI008:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI010:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI012:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
str(spatial.dat[[1]]@MASKS, 3)
## List of 1
## $ cell.mask:List of 1
## ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
###################################################################################
### Rename rasters (if required)
###################################################################################
Here you can check the consistency of the raster names, and if required, adjust those names.
### Check channel names
channel.names <- list()
for(i in names(spatial.dat)){
channel.names[[i]] <- names(spatial.dat[[i]]@RASTERS)
}
t(as.data.frame(channel.names))
## [,1] [,2] [,3] [,4] [,5]
## ROI002 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI004 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI006 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI008 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI010 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI012 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## [,6] [,7] [,8]
## ROI002 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI004 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI006 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI008 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI010 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI012 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
### List of corrections (first entry is the 'correct' one)
# corrections <- list(c('CD4','Cd4'),
# c('CD8','CD8a')
# )
### Replace the 'incorrect' names
# for(i in names(spatial.dat)){
# # i <- names(spatial.dat)[[1]]
#
# for(a in c(1:length(corrections))){
# # a <- 1
#
# trg <- which(names(spatial.dat[[i]]@RASTERS) == corrections[[a]][2])
# if(length(trg) != 0){
# names(spatial.dat[[i]]@RASTERS)[trg] <- corrections[[a]][1]
# }
# }
# }
### Check channel names
# channel.names <- list()
#
# for(i in names(spatial.dat)){
# channel.names[[i]] <- names(spatial.dat[[i]]@RASTERS)
# }
#
# t(as.data.frame(channel.names))
###################################################################################
### Generate polygons and outlines
###################################################################################
### Generate polygons and outlines
for(i in names(mask.types)){
spatial.dat <- do.create.outlines(dat = spatial.dat, mask.name = i)
}
### Checks
str(spatial.dat, 3)
## List of 6
## $ ROI002:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI004:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI006:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI008:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI010:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
## $ ROI012:Formal class 'spatial' [package "Spectre"] with 3 slots
## .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
## .. ..@ MASKS :List of 1
## .. ..@ DATA : list()
str(spatial.dat[[1]]@MASKS, 2)
## List of 1
## $ cell.mask:List of 4
## ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
## ..$ polygons :Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..$ outlines :'data.frame': 93518 obs. of 7 variables:
## ..$ centroids :Formal class 'SpatialPoints' [package "sp"] with 3 slots
###################################################################################
### Mask QC plots
###################################################################################
Here you can choose a ‘base’ raster to use for the spatial plot, and which mask you would like to plot.
### Mask plot setup
setwd(OutputDirectory)
dir.create('Plots - cell masks')
setwd('Plots - cell masks')
as.matrix(names(spatial.dat[[1]]@RASTERS))
## [,1]
## [1,] "CD11b_Sm149"
## [2,] "CD20_Dy161"
## [3,] "CD3_Er170"
## [4,] "CD4_Gd156"
## [5,] "CD45_Sm152"
## [6,] "CD8a_Dy162"
## [7,] "DNA1_Ir191"
## [8,] "DNA3_Ir193"
base <- 'DNA1_Ir191'
base
## [1] "DNA1_Ir191"
as.matrix(names(spatial.dat[[1]]@MASKS))
## [,1]
## [1,] "cell.mask"
mask <- 'cell.mask'
mask
## [1] "cell.mask"
### Create plots
for(i in names(spatial.dat)){
make.spatial.plot(dat = spatial.dat,
image.roi = i,
image.channel = base,
mask.outlines = mask)
}
You will see plots that look like this for each ROI, allowing you to assess the suitability of the mask.
### Create plots
for(i in names(spatial.dat)[1]){
make.spatial.plot(dat = spatial.dat,
image.roi = i,
image.channel = base,
mask.outlines = mask)
}
###################################################################################
### Calculate cellular data and plot
###################################################################################
### Calculate cellular data for each cell mask (this step may take some time)
spatial.dat <- do.extract(spatial.dat, 'cell.mask')
str(spatial.dat, 3)
spatial.dat[[1]]@DATA
### Factor plot setup
setwd(OutputDirectory)
dir.create('Plots - factors')
setwd('Plots - factors')
as.matrix(names(spatial.dat))
## [,1]
## [1,] "ROI002"
## [2,] "ROI004"
## [3,] "ROI006"
## [4,] "ROI008"
## [5,] "ROI010"
## [6,] "ROI012"
plot.rois <- names(spatial.dat)[c(1:2)]
plot.rois
## [1] "ROI002" "ROI004"
as.matrix(names(spatial.dat[[1]]@RASTERS))
## [,1]
## [1,] "CD11b_Sm149"
## [2,] "CD20_Dy161"
## [3,] "CD3_Er170"
## [4,] "CD4_Gd156"
## [5,] "CD45_Sm152"
## [6,] "CD8a_Dy162"
## [7,] "DNA1_Ir191"
## [8,] "DNA3_Ir193"
base <- 'DNA1_Ir191'
base
## [1] "DNA1_Ir191"
plot.factors <- names(spatial.dat[[1]]@MASKS)[-which('cell.mask' == names(spatial.dat[[1]]@MASKS))]
plot.factors
## character(0)
plot.exp <- names(spatial.dat[[1]]@RASTERS)
plot.exp
## [1] "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## [6] "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
### Make factor plots
for(i in plot.rois){
setwd(OutputDirectory)
setwd('Plots - factors')
dir.create(i)
setwd(i)
for(a in plot.factors){
make.spatial.plot(dat = spatial.dat,
image.roi = i,
image.channel = base,
mask.outlines = mask,
cell.dat = 'CellData',
cell.col = a,
cell.col.type = 'factor',
title = paste0(a, ' - ', i))
}
}
### Make exp plots
for(i in plot.rois){
setwd(OutputDirectory)
setwd('Plots - factors')
dir.create(i)
setwd(i)
for(a in plot.exp){
make.spatial.plot(dat = spatial.dat,
image.roi = i,
image.channel = base,
mask.outlines = mask,
cell.dat = 'CellData',
cell.col = a,
title = paste0(a, ' - ', i))
}
}
###################################################################################
### Extract cellular data and annotate
###################################################################################
### Extract cellular data
cell.dat <- do.pull.data(spatial.dat, 'CellData')
cell.dat
## ROI ID x y Area CD11b_Sm149 CD20_Dy161 CD3_Er170
## 1: ROI002 1 276.24140 218.403834 100982 0.3623715 6.8142638 0.4266008
## 2: ROI002 2 222.54000 3.200000 50 0.2800000 1.9000000 0.4000000
## 3: ROI002 3 247.44643 3.660714 56 0.3928571 0.4464286 0.1428571
## 4: ROI002 4 333.27941 4.250000 68 0.3529412 1.2205882 0.4558823
## 5: ROI002 5 388.37755 4.785714 49 0.1428571 0.5714286 0.1632653
## ---
## 15564: ROI012 1423 216.41892 494.986486 37 0.2972973 0.2972973 0.2432432
## 15565: ROI012 1424 59.82099 495.919753 81 0.3333333 1.0000000 0.2592593
## 15566: ROI012 1425 432.20732 495.731707 82 0.6829268 0.8780488 0.2560976
## 15567: ROI012 1426 397.91463 495.548780 41 0.3170732 0.5121951 0.4146341
## 15568: ROI012 1427 420.92105 496.578947 38 0.2105263 1.1578947 2.6842105
## CD4_Gd156 CD45_Sm152 CD8a_Dy162 DNA1_Ir191 DNA3_Ir193
## 1: 0.6018696 3.1577609 1.5598127 42.90354 79.08650
## 2: 0.6000000 1.5800000 4.0999999 44.28000 84.58000
## 3: 0.2678571 0.6607143 3.0357144 29.41072 55.87500
## 4: 1.0441177 2.3529413 3.3823528 29.25000 54.79412
## 5: 0.1224490 0.3265306 2.4489796 43.75510 81.75510
## ---
## 15564: 0.6216216 1.0270270 0.3243243 40.29730 76.18919
## 15565: 0.7283950 0.6790124 1.8148148 29.44444 55.50617
## 15566: 1.0487804 2.0975609 0.9878049 42.02439 77.53658
## 15567: 0.9756098 3.1219511 1.2439024 34.85366 63.07317
## 15568: 0.9210526 10.2368422 4.3421054 33.28947 60.76316
###################################################################################
### Area calculation
###################################################################################
area.totals <- do.calculate.area(spatial.dat)
area.totals
## ROI Total
## 1: ROI002 500.4998
## 2: ROI004 500.4998
## 3: ROI006 500.0000
## 4: ROI008 500.4998
## 5: ROI010 500.4998
## 6: ROI012 500.4998
###################################################################################
### Save data
###################################################################################
### Output QS and CSV file
setwd(OutputDirectory)
dir.create('Data')
setwd('Data')
qsave(spatial.dat, "spatial.dat.qs")
fwrite(cell.dat, 'cell.dat.csv')
fwrite(area.totals, 'area.totals.csv')
### Pull cellular data and write FCS file from each ROI independently
setwd(OutputDirectory)
dir.create('FCS files')
setwd('FCS files')
for(i in names(spatial.dat)){
## Extract data and setup cols
tmp <- list()
tmp[[i]] <- spatial.dat[[i]]
cell.dat <- do.pull.data(tmp, 'CellData')
cell.dat <- do.asinh(cell.dat, names(spatial.dat[[i]]@RASTERS), cofactor = 1)
### Invert y axis
all.neg <- function(test) -1*abs(test)
y_invert <- cell.dat[['y']]
y_invert <- all.neg(y_invert)
cell.dat[['y_invert']] <- y_invert
### Write FCS files
write.files(cell.dat, i, write.csv = FALSE, write.fcs = TRUE)
rm(cell.dat)
rm(i)
}
Coming soon!
Coming soon!